1 Overview

This project aims to connect floral phenology to reproductive success on the individual level for three subalpine plants. Data were collected at the Rocky Mountain Biological Laboratory during the 2019 field season. Focal species include Mertensia fusiformis, Potentilla pulcherrima, and Delphinium nuttallianum.

Our questions specifically ask: 1. how does blooming time affect seed set for individiuals in a population? 2. how do abiotic factors such as soil moisture affect individual seed set? 3. how do biotic factors such as conspecific density and pollen limitation affect individual seed set?

We tagged individuals of the three species in control and manipulated plots in eight sites across East River valley and Washington Gulch. We tracked blooming time of those tagged individuals and conducted a pollen limitation experiment with two treatments, open and pollen-supplemented. Seed set was quantified as the count of total developed seeds and proportion of developed seeds in each individual.

Because there were sources of non-independence in this study, our statistical approach involves generalized linear mixed effects models. Our models includes site as a random effect and the following fixed effects: the interaction between being before or after the population peak and deviation from the peak and the interaction between conspecific density and plot treatment. A separate model included pollen limitation treatment.

2 Load Packages

rm(list = ls()) #clearing environment
suppressWarnings(suppressPackageStartupMessages(require(knitr)))
suppressWarnings(suppressPackageStartupMessages(require(tidyverse)))
suppressWarnings(suppressPackageStartupMessages(require(lme4)))
suppressWarnings(suppressPackageStartupMessages(require(ggplot2)))
suppressWarnings(suppressPackageStartupMessages(require(RColorBrewer)))
suppressWarnings(suppressPackageStartupMessages(require(fitdistrplus)))
suppressWarnings(suppressPackageStartupMessages(require(glmmADMB)))
suppressWarnings(suppressPackageStartupMessages(require(MuMIn)))
suppressWarnings(suppressPackageStartupMessages(require(lubridate)))
suppressWarnings(suppressPackageStartupMessages(require(DHARMa)))
suppressWarnings(suppressPackageStartupMessages(require(glmmTMB)))
suppressWarnings(suppressPackageStartupMessages(require(dplyr)))
suppressWarnings(suppressPackageStartupMessages(require(stringr)))
suppressWarnings(suppressPackageStartupMessages(require(scales)))
suppressWarnings(suppressPackageStartupMessages(require(lmerTest)))

3 Data Import & Basic Formatting

The imported csv files contain the following for each individual of a species: bloom start and end estimates, relative position of blooming in the community, soil moisture variables, and seed counts. These elements were compiled from five different data sets and calculated in a separate document. To view the data cleaning process, refer to the data cleaning R Markdown file, Data_Cleaning.Rmd.

3.1 Mertensia Data

The first dataset contains the individual data for Mertensia fusiformis, which was found in four of our sites. For Mertensia, we calculated the proportion of developed seeds in addition to the seed counts. We also reformatted the relative position as a deviation variable (the absolute value of the relative position) and a early/late variable (the direction of the relative position).

mert = read.csv("Mertensia_final2.csv") # read in Mertensia data
mert$prop.dev.seeds<-(mert$dev.seed/mert$total.seeds) #adding proportion of developed seeds
mert$deviation<-abs(mert$relative.position) #getting deviation from the peak

mert$early.late<-factor(NA, levels = c("early","late")) #creating early/late variables based on relative position
mert$early.late[mert$relative.position<0]<-"early"
mert$early.late[mert$relative.position>0]<-"late"
mert<-mert[!is.na(mert$early.late),]

mert<-mert[!is.na(mert$dev.seed),]
mert$prop.dev.seeds[is.na(mert$prop.dev.seeds)] <-0

3.2 Delphinium Data

The second dataset contains individual data for Delphinium nuttallianum. Delphinium was found in four of our sites. For Delphinium, we calculated proportion, deviation, and early/late value like we did with the Mertensia data.

delph = read.csv("Delphinium_final2.csv") # read in Delphinium data
delph$undev.seed1[is.na(delph$undev.seed1)] <- 0 #converting NA to 0
delph$dev.seed.total<-delph$dev.seed+delph$dev.seed1 #creating developed seeds column
delph$undev.seed.total<-delph$undev.seed+delph$undev.seed1 #creating undeveloped seeds column
delph$prop.dev.seeds<-(delph$dev.seed.total/(delph$dev.seed.total+delph$undev.seed.total)) #adding proportion of developed seeds

delph$deviation<-abs(delph$relative.position) #getting deviation from the peak
delph$early.late<-factor(NA, levels = c("early","late")) #creating early/late variables based on relative position
delph$early.late[delph$relative.position<0]<-"early"
delph$early.late[delph$relative.position>0]<-"late"
delph<-delph[!is.na(delph$early.late),]

delph<-delph[!is.na(delph$dev.seed.total),]
delph$prop.dev.seeds[is.na(delph$prop.dev.seeds)] <-0

3.3 Potentilla Data

The third dataset contains Potentilla pulcherrima data. Potentilla was set up in all but one of our sites. Potentilla data only includes total seed counts because we couldn't count the undeveloped seeds. However, we calculated the deviation and early/late values like the previous species.

pot = read.csv("Potentilla_final2.csv") # read in Potentilla data
pot$deviation<-abs(pot$relative.position) #getting deviation from the peak

pot$early.late<-factor(NA, levels = c("early","late")) #creating early/late variables based on relative position
pot$early.late[pot$relative.position<0]<-"early"
pot$early.late[pot$relative.position>0]<-"late"
pot<-pot[!is.na(pot$early.late),] #removing NA values in early/late column

Because many of the Potentilla individuals were still in bloom when we finished our field season, we had to calculate a midpoint and relative position based on those individuals that had senesced. A key factor that determines duration of blooming is the number of flowers on the individual. The data imported includes those relative position estimates for the individuals that didn't senesce before we finished our field season. For more information on how those numbers were calculated, see Data_Cleaning.Rmd.

3.4 Flower Count Data

The community flower count data is essential to the number of conspecifics effect. Each of the species has to be analyzed separately.

flower.counts=read.csv("flower_counts.csv") #read in flower count data

3.4.1 Mertensia

We calculated the average number of Mertensia flowers in each site and plot in each week. This number was then merged with the Mertensia individuals to add the number of conspecifics variable.

mert.flower.counts<-flower.counts[(flower.counts$plant=="Mertensia fusiformis"),] #limit flower counts to Mertensia data
mert.counts.by.site<-aggregate(x=mert.flower.counts$total_flowers,by=list(mert.flower.counts$date,mert.flower.counts$site,mert.flower.counts$plot),FUN="mean") #calculate mean flower counts by week and site
colnames(mert.counts.by.site)<-c("date","site","plot.treat","number.conspecifics")

mert.counts.by.site$site<-as.character(mert.counts.by.site$site) #changing class to adjust names
mert.counts.by.site$plot.treat<-as.character(mert.counts.by.site$plot.treat)

mert.counts.by.site$site[mert.counts.by.site$site=="Avery Picnic"]<-"AveryPicnic" #fixing names
mert.counts.by.site$site[mert.counts.by.site$site=="Stupid Falls"]<-"StupidFalls"
mert.counts.by.site$site[mert.counts.by.site$site=="Wash 3C"]<-"WG3C"
mert.counts.by.site$plot.treat[!(mert.counts.by.site$plot.treat=="Kaysee")]
mert.counts.by.site$plot.treat[mert.counts.by.site$plot.treat=="Control"]<-"control"
mert.counts.by.site$plot.treat[mert.counts.by.site$plot.treat=="Manipulated"]<-"manip"

mert$week.estimate<-strftime(mert$midpoint, format = "%V") #calculating week numbers in preparation for the merge
mert.counts.by.site$week.estimate<-strftime(mert.counts.by.site$date, format = "%V")

mert<-merge(mert,mert.counts.by.site, by.x = c("site","plot.treat","week.estimate"), by.y = c("site","plot.treat","week.estimate")) #merging by site, plot treatment, and week estimate

3.4.2 Delphinium

The same was repeated for Delphinium.

delph.flower.counts<-flower.counts[(flower.counts$plant=="Delphinium nuttallianum"),] #limit flower counts to Delphinium data
delph.counts.by.site<-aggregate(x=delph.flower.counts$total_flowers,by=list(delph.flower.counts$date,delph.flower.counts$site,delph.flower.counts$plot),FUN="mean") #calculate mean flower counts by week and site
colnames(delph.counts.by.site)<-c("date","site","plot.treat","number.conspecifics")

delph.counts.by.site$site<-as.character(delph.counts.by.site$site) #changing class to adjust names
delph.counts.by.site$plot.treat<-as.character(delph.counts.by.site$plot.treat)

delph.counts.by.site$site[delph.counts.by.site$site=="Rustlers Gulch"]<-"RustlersGulch" #fixing names
delph.counts.by.site$site[delph.counts.by.site$site=="Wash 3C"]<-"WG3C"
#delph.counts.by.site$plot.treat[!(delph.counts.by.site$plot.treat=="Kaysee")]
delph.counts.by.site$plot.treat[delph.counts.by.site$plot.treat=="Control"]<-"control"
delph.counts.by.site$plot.treat[delph.counts.by.site$plot.treat=="Manipulated"]<-"manip"

delph$week.estimate<-strftime(delph$midpoint, format = "%V") #calculating week numbers in preparation for the merge
delph.counts.by.site$week.estimate<-strftime(delph.counts.by.site$date, format = "%V")

delph<-merge(delph,delph.counts.by.site, by.x = c("site","plot.treat","week.estimate"), by.y = c("site","plot.treat","week.estimate")) #merging by site, plot treatment, and week estimate

3.4.3 Potentilla

The same was done for Potentilla.

pot.flower.counts<-flower.counts[(flower.counts$plant=="Potentilla pulcherrima"),] #limit flower counts to Potentilla data
pot.counts.by.site<-aggregate(x=pot.flower.counts$total_flowers,by=list(pot.flower.counts$date,pot.flower.counts$site,pot.flower.counts$plot),FUN="mean") #calculate mean flower counts by week and site
colnames(pot.counts.by.site)<-c("date","site","plot.treat","number.conspecifics")

pot.counts.by.site$site<-as.character(pot.counts.by.site$site) #changing class to adjust names
pot.counts.by.site$plot.treat<-as.character(pot.counts.by.site$plot.treat)

pot.counts.by.site$site[pot.counts.by.site$site=="Rustlers Gulch"]<-"RustlersGulch" #fixing names
pot.counts.by.site$site[pot.counts.by.site$site=="Wash 3C"]<-"WG3C"
pot.counts.by.site$site[pot.counts.by.site$site=="Bellview Bench"]<-"BellviewBench"
pot.counts.by.site$site[pot.counts.by.site$site=="Avery Picnic"]<-"AveryPicnic"
pot.counts.by.site$site[pot.counts.by.site$site=="Stupid Falls"]<-"StupidFalls"
#pot.counts.by.site$plot.treat[!(pot.counts.by.site$plot.treat=="Kaysee")]
pot.counts.by.site$plot.treat[pot.counts.by.site$plot.treat=="Control"]<-"control"
pot.counts.by.site$plot.treat[pot.counts.by.site$plot.treat=="Manipulated"]<-"manip"

pot$week.estimate<-strftime(pot$midpoint, format = "%V") #calculating week numbers in preparation for the merge
pot.counts.by.site$week.estimate<-strftime(pot.counts.by.site$date, format = "%V")

pot<-merge(pot,pot.counts.by.site, by.x = c("site","plot.treat","week.estimate"), by.y = c("site","plot.treat","week.estimate")) #merging by site, plot treatment, and week estimate

4 Mixed Effects Models with Count Data

Generalized linear mixed effects models were appropriate for this analysis because the data contain sources of non-independence. Individual plants were located in the same site as other individuals and were likely genetically similar. Therefore, we considered site a random effect.

The response variable is total seeds, which is discrete and non-zero. The interaction between deviation and the early/late factor was a fixed effect, as well as the interaction between plot treatment and conspecific density. We would have included soil moisture as a covariate and fixed effect in our model, but soil moisture was not important for seed production. Each species is analyzed separately.

4.1 Mertensia

We started by limiting the individuals to the open treatment (excluding the hand-pollinated treatment) and by checking the distribution of the total seed counts.

mert.open<-mert[(mert$treat=="open"),] #removed hand pollination treatment
mert.nbinom<-fitdist(mert.open$dev.seed, "nbinom") #fitting Mertensia data to negative binomial distribution using fitdistplus
plot(mert.nbinom) #plotting to see the fit

The Mertensia seed data fits a negative binomial distribution. The glmmTMB package is a good package for modeling data with a negative binomial. The glmmTMB package is also suitable for data that could be overdispersed or zero-inflated.

4.1.1 Main Model

#mert.glmmtmb<-glmmTMB(dev.seed ~ flowers:(deviation*early.late) + flowers:number.conspecifics/plot.treat  + (1|site), family = "nbinom2", data = mert.open) #putting count data into glmmTMB for DHARMa package
#summary(mert.glmmtmb) #summary of glmmTMB model

#mert.counts.simulation<-simulateResiduals(fittedModel = mert.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
#plot(mert.counts.simulation) #plotting simulation
#testDispersion(mert.counts.simulation) #checking for overdispersion
#testZeroInflation(mert.counts.simulation) #checking for zero inflation

We detected a significant effect of early vs. late blooming on the total number of seeds in Mertensia individuals, but when we checked for overdispersion and zero-inflation in the data using the DHARMa package, the data were slightly zero-inflated. To address this, we adjusted our glmmtmb model to include zero-inflation.

mert.glmmtmb<-glmmTMB(dev.seed ~ flowers:(deviation*early.late) + flowers:number.conspecifics/plot.treat  + (1|site), ziformula=~1, family = "nbinom2", data = mert.open) #putting count data into glmmTMB for DHARMa package
summary(mert.glmmtmb) #summary of glmmTMB model
##  Family: nbinom2  ( log )
## Formula:          
## dev.seed ~ flowers:(deviation * early.late) + flowers:number.conspecifics/plot.treat +  
##     (1 | site)
## Zero inflation:            ~1
## Data: mert.open
## 
##      AIC      BIC   logLik deviance df.resid 
##    386.2    402.4   -184.1    368.2       36 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.07728  0.278   
## Number of obs: 45, groups:  site, 4
## 
## Overdispersion parameter for nbinom2 family (): 3.03 
## 
## Conditional model:
##                                               Estimate Std. Error z value
## (Intercept)                                  2.613e+00  3.761e-01   6.949
## flowers:deviation                            2.433e-03  1.540e-03   1.579
## flowers:early.latelate                       1.931e-02  1.009e-02   1.914
## flowers:number.conspecifics                 -2.368e-05  8.772e-05  -0.270
## flowers:deviation:early.latelate            -2.693e-03  2.934e-03  -0.918
## flowers:number.conspecifics:plot.treatmanip  2.973e-04  1.937e-04   1.535
##                                             Pr(>|z|)    
## (Intercept)                                 3.68e-12 ***
## flowers:deviation                             0.1143    
## flowers:early.latelate                        0.0557 .  
## flowers:number.conspecifics                   0.7872    
## flowers:deviation:early.latelate              0.3587    
## flowers:number.conspecifics:plot.treatmanip   0.1248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.3800     0.5513  -4.317 1.58e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mert.counts.simulation<-simulateResiduals(fittedModel = mert.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.counts.simulation) #plotting simulation

testDispersion(mert.counts.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.90279, p-value = 0.768
## alternative hypothesis: two.sided
testZeroInflation(mert.counts.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0417, p-value = 1
## alternative hypothesis: two.sided

When number of flowers is taken into account, seed set for Mertensia individuals may be affected by blooming before or after the population peak. Zero-inflation is accounted for in the new model, and the model for Mertensia count data is no longer zero-inflated.

To figure out the direction of the effect, we looked at the plots (shown below). Blooming after the population peak may increase seed set for Mertensia individuals.

4.1.2 Pollen Limitation Model

To test the pollen limitation individuals, we created a separate fixed effects model. This model has site and plot treatment as random effects and the interaction between pollen limitation treatment and number of flowers as a fixed effect. We included this interaction because the difference in seed set between the open and hand treatments is heavily dependent on the number of flowers. The interaction accounts for the wide variation of number of flowers and ultimately developed seeds. The response variable is still total number of developed seeds in Mertensia individuals.

#mert.treat.glmmtmb<-glmmTMB(dev.seed ~ treat/flowers + (1|site/plot.treat), family = "nbinom2", data = mert) #Mertensia model for negative binomial data with glmmTMB package
#summary(mert.treat.glmmtmb) #summary of model output

#mert.counts.treat.simulation<-simulateResiduals(fittedModel = mert.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
#plot(mert.counts.treat.simulation) #plotting simulation
#testDispersion(mert.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(mert.counts.treat.simulation) #checking for zero inflation

We detected an effect of pollen limitation treatment on Mertensia total developed seeds, but we had to test the model assumptions of our pollen limitation model as well. The data in the pollen limitation model for Mertensia are not overdispersed, but they are zero-inflated. We then adjusted our model to include zero-inflation and ran the overdispersion and zero-inflation tests again.

mert.treat.glmmtmb<-glmmTMB(dev.seed ~ treat/flowers + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = mert) #Mertensia model for negative binomial data with glmmTMB package
summary(mert.treat.glmmtmb) #summary of model output
##  Family: nbinom2  ( log )
## Formula:          dev.seed ~ treat/flowers + (1 | site/plot.treat)
## Zero inflation:            ~1
## Data: mert
## 
##      AIC      BIC   logLik deviance df.resid 
##    750.6    770.9   -367.3    734.6       85 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance Std.Dev.
##  plot.treat:site (Intercept) 0.03016  0.1737  
##  site            (Intercept) 0.07045  0.2654  
## Number of obs: 93, groups:  plot.treat:site, 8; site, 4
## 
## Overdispersion parameter for nbinom2 family (): 5.64 
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            2.372696   0.244703   9.696  < 2e-16 ***
## treatopen-hand         0.011735   0.267892   0.044    0.965    
## treatopen:flowers      0.024131   0.004452   5.420 5.95e-08 ***
## treatopen-hand:flowers 0.023539   0.004557   5.165 2.40e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.8836     0.4667  -6.178 6.49e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mert.counts.treat.simulation<-simulateResiduals(fittedModel = mert.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.counts.treat.simulation) #plotting simulation

testDispersion(mert.counts.treat.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.87794, p-value = 0.648
## alternative hypothesis: two.sided
testZeroInflation(mert.counts.treat.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0238, p-value = 1
## alternative hypothesis: two.sided

The data are no longer zero-inflated, and we still detected an effect of pollen limitation on the Mertensia seed counts. Mertensia individuals are pollen limited.

4.1.3 Figures

Below is the plot of seed set vs. relative position of blooming in the population for Mertensia individuals. Plot treatment is indicated by individual color.

ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$dev.seed, group=mert.open$plot.treat)) +
  labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Blooming Time on Total Developed Seeds") +
  geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
  geom_smooth()

ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$dev.seed, group=mert.open$early.late)) +
  labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Blooming Time on Total Developed Seeds") +
  theme_classic() +
  geom_point(aes(color=early.late),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
  geom_smooth(method = "glm", aes(color=early.late))

We then plotted the interaction between deviation and early/late instead of relative position in order to fit a linear model.

ggplot(mert.open, aes(x=mert.open$deviation, y=mert.open$dev.seed/mert.open$flowers, group= mert.open$early.late)) +
  labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Deviation from Peak on Total Developed Seeds") + theme_classic() +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
  geom_smooth(method="glm",aes(color=early.late)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))

Then we plotted the effect of conspecific density and plot treatment on total developed seeds.

ggplot(mert.open, aes(x=mert.open$number.conspecifics, y=mert.open$dev.seed/mert.open$flowers, group= mert.open$plot.treat)) +
  labs(y="Total Developed Seeds per Flower",x="Conspecific Density", title="Effect of Mertensia Conspecific Density on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
  geom_smooth(method="glm",aes(color=plot.treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))

Finally, we plotted all of the data with the pollinator active period. The pollinator active period was the time during which we saw interactions between pollinators and our species. The active period is shaded.

mert$midpoint2<-format(as.Date(mert$midpoint), format="%m-%d")

ggplot(mert, aes(x=mert$midpoint2, y=mert$dev.seed/mert$flowers, group= mert$treat)) +
  labs(y="Total Developed Seeds per Flower",x="Date", title="Effect of Mertensia Blooming Time and Pollinator Activity on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=treat)) +
  geom_smooth(aes(color=treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Open or Hand-pollinated") +
  scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  annotate("rect",xmin="06-19", xmax="07-03",ymin=0,ymax=Inf,alpha=0.1)

4.2 Delphinium

For Delphinium, we again limited the individuals to open treatment individuals.

delph.open<-delph[(delph$treat=="open"),] #removed hand pollination treatment
delph.nbinom<-fitdist(delph.open$total.seeds, "nbinom") #fitting negative binomial distribution using fitdistplus package
plot(delph.nbinom) #plotting distribution fit

The Delphinium data fits a negative binomial distribution. We used the glmmTMB package again.

4.2.1 Main Model

delph.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation*early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = delph.open) #Delphinium model for negative binomial data with glmmTMB package
#removed treat fixed effect (only open)
summary(delph.glmmtmb) #summary of model output
##  Family: nbinom2  ( log )
## Formula:          
## total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat +  
##     (1 | site)
## Data: delph.open
## 
##      AIC      BIC   logLik deviance df.resid 
##    403.9    417.6   -193.9    387.9       33 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.02556  0.1599  
## Number of obs: 41, groups:  site, 4
## 
## Overdispersion parameter for nbinom2 family (): 1.37 
## 
## Conditional model:
##                                                Estimate Std. Error z value
## (Intercept)                                    2.990040   0.369707   8.088
## X.flowers:deviation                            0.026479   0.011771   2.249
## X.flowers:early.latelate                       0.137212   0.168217   0.816
## X.flowers:number.conspecifics                  0.001297   0.001946   0.666
## X.flowers:deviation:early.latelate            -0.010662   0.036032  -0.296
## X.flowers:number.conspecifics:plot.treatmanip -0.002071   0.001801  -1.150
##                                               Pr(>|z|)    
## (Intercept)                                   6.09e-16 ***
## X.flowers:deviation                             0.0245 *  
## X.flowers:early.latelate                        0.4147    
## X.flowers:number.conspecifics                   0.5051    
## X.flowers:deviation:early.latelate              0.7673    
## X.flowers:number.conspecifics:plot.treatmanip   0.2501    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
delph.counts.simulation<-simulateResiduals(fittedModel = delph.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
plot(delph.counts.simulation) #plotting simulation

testDispersion(delph.counts.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.79078, p-value = 0.448
## alternative hypothesis: two.sided
testZeroInflation(delph.counts.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 4.2735, p-value = 0.136
## alternative hypothesis: two.sided

When flower number is taken into account, deviation from the population peak has a significant effect on Delphinium developed seed count. We tested the model assumptions as we did with Mertensia seed count data. The model for the Delphinium count data is not overdispersed or zero-inflated.

To determine the direction of deviation effect, we looked at the plot. When flower number is taken into account, seed set of Delphinium individuals increase with deviation from the population peak.

4.2.2 Pollen Limitation Model

Again, we tested the pollen limitation in a separate mixed effects model.

#delph.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), family = "nbinom2", data = delph) #Mertensia model for negative binomial data with glmmTMB package
#summary(delph.treat.glmmtmb) #summary of model output

#delph.counts.treat.simulation<-simulateResiduals(fittedModel = delph.treat.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
#plot(delph.counts.treat.simulation) #plotting simulation
#testDispersion(delph.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(delph.counts.treat.simulation) #checking for zero inflation

We detected an effect of pollen limitation on the total developed seeds of Delphinium individuals. We again checked for overdispersion and zero-inflation in the pollen limitation data. The Delphinum pollen limitation model is zero-inflated. We adjusted our model to include the zero-inflation formula.

delph.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = delph) #Mertensia model for negative binomial data with glmmTMB package
summary(delph.treat.glmmtmb) #summary of model output
##  Family: nbinom2  ( log )
## Formula:          total.seeds ~ treat/X.flowers + (1 | site/plot.treat)
## Zero inflation:               ~1
## Data: delph
## 
##      AIC      BIC   logLik deviance df.resid 
##    811.2    830.6   -397.6    795.2       76 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance  Std.Dev. 
##  plot.treat:site (Intercept) 2.415e-09 4.914e-05
##  site            (Intercept) 6.478e-02 2.545e-01
## Number of obs: 84, groups:  plot.treat:site, 7; site, 4
## 
## Overdispersion parameter for nbinom2 family (): 1.97 
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               2.51490    0.40312   6.239 4.42e-10 ***
## treatopen-hand            0.69778    0.50148   1.391 0.164096    
## treatopen:X.flowers       0.27607    0.07661   3.603 0.000314 ***
## treatopen-hand:X.flowers  0.16104    0.06725   2.395 0.016636 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.1192     0.5794  -5.383 7.31e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
delph.counts.treat.simulation<-simulateResiduals(fittedModel = delph.treat.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
plot(delph.counts.treat.simulation) #plotting simulation

testDispersion(delph.counts.treat.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.79907, p-value = 0.352
## alternative hypothesis: two.sided
testZeroInflation(delph.counts.treat.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.065, p-value = 1
## alternative hypothesis: two.sided

The new model is not zero-inflated. We still see an effect of pollen limitation on the Delphinium count data. Delphinium individuals are pollen limited.

4.2.3 Figures

Below are the plots for the Delphinium individuals.

ggplot(delph.open, aes(x=delph.open$relative.position, y=delph.open$total.seeds, group=plot.treat)) +
  labs(title="Effect of Delphinium Blooming Time on Total Developed Seeds", y="Total Developed Seeds",x="Days Relative to Population Peak Bloom") +
  geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
  geom_smooth()

ggplot(delph.open, aes(x=delph.open$relative.position, y=delph.open$total.seeds, group=delph.open$early.late)) +
  labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Delphinium Blooming Time on Total Developed Seeds") +
  theme_classic() +
  geom_point(aes(color=early.late),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
  geom_smooth(method = "glm", aes(color=early.late))

ggplot(delph.open, aes(x=delph.open$deviation, y=(delph.open$total.seeds/delph.open$X.flowers), group= delph.open$early.late)) +
  labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom", title="Effect of Delphinium Deviation from Peak on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
  geom_smooth(method="glm",aes(color=early.late)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))

ggplot(delph.open, aes(x=delph.open$number.conspecifics, y=delph.open$total.seeds/delph.open$X.flowers, group= delph.open$plot.treat)) +
  labs(y="Total Developed Seeds per Flower",x="Conspecific Density", title="Effect of Delphinium Conspecific Density on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
  geom_smooth(method="glm",aes(color=plot.treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))

Finally, we plotted the pollinator active period with the total number of developed seeds and midpoint of individuals.

delph$midpoint2<-format(as.Date(delph$midpoint), format="%m-%d")

ggplot(delph, aes(x=delph$midpoint2, y=delph$total.seeds/delph$X.flowers, group= delph$treat)) +
  labs(y="Total Developed Seeds",x="Date", title="Effect of Delphinium Blooming and Pollinator Activity on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=treat)) +
  geom_smooth(aes(color=treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Open or Hand-pollinated") +
  scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  annotate("rect",xmin="06-27", xmax="07-18",ymin=0,ymax=Inf,alpha=0.1)

4.3 Potentilla

We began by limiting the data to the open treatment Potentilla individuals.

pot.open<-pot[(pot$treat=="open"),] #removed hand pollination treatment
pot.nbinom<-fitdist(pot.open$total.seeds, "nbinom") #fitting negative binomial distribution using fitdistplus package
plot(pot.nbinom) #plotting distribution fit

The Potentilla seed data fits a negative binomial distribution. We used the glmmTMB package again.

4.3.1 Main Model

#pot.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = pot.open) #Potentilla model for negative binomial data with glmmTMB package
#summary(pot.glmmtmb) #summary of model output

#pot.counts.simulation<-simulateResiduals(fittedModel = pot.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
#plot(pot.counts.simulation) #plotting simulation
#testDispersion(pot.counts.simulation) #checking for overdispersion
#testZeroInflation(pot.counts.simulation) #checking for zero inflation

The Potentilla data are not overdispersed, but they are zero-inflated. To fix this, we had to include zero inflation in our model and again check for zero-inflation.

pot.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), ziformula=~1, family = "nbinom2", data = pot.open) #new Potentilla model for negative binomial data with glmmTMB package
summary(pot.glmmtmb) #summary of model output
##  Family: nbinom2  ( log )
## Formula:          
## total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat +  
##     (1 | site)
## Zero inflation:               ~1
## Data: pot.open
## 
##      AIC      BIC   logLik deviance df.resid 
##    280.7    291.4   -131.4    262.7       15 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 3.537e-10 1.881e-05
## Number of obs: 24, groups:  site, 4
## 
## Overdispersion parameter for nbinom2 family (): 2.95 
## 
## Conditional model:
##                                                 Estimate Std. Error z value
## (Intercept)                                    3.657e+00  4.095e-01   8.930
## X.flowers:deviation                            1.220e-02  4.421e-03   2.760
## X.flowers:early.latelate                       7.605e-02  2.733e-02   2.783
## X.flowers:number.conspecifics                  1.631e-05  3.236e-04   0.050
## X.flowers:deviation:early.latelate            -1.423e-02  4.629e-03  -3.074
## X.flowers:number.conspecifics:plot.treatmanip  4.595e-04  3.227e-04   1.424
##                                               Pr(>|z|)    
## (Intercept)                                    < 2e-16 ***
## X.flowers:deviation                            0.00578 ** 
## X.flowers:early.latelate                       0.00539 ** 
## X.flowers:number.conspecifics                  0.95981    
## X.flowers:deviation:early.latelate             0.00211 ** 
## X.flowers:number.conspecifics:plot.treatmanip  0.15444    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   -3.138      1.024  -3.065  0.00218 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pot.counts.simulation<-simulateResiduals(fittedModel = pot.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
plot(pot.counts.simulation) #plotting simulation

testDispersion(pot.counts.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.8047, p-value = 0.544
## alternative hypothesis: two.sided
testZeroInflation(pot.counts.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0776, p-value = 1
## alternative hypothesis: two.sided

When we incorporated zero-inflation into the model, we found a significant effect of deviation from the population peak, blooming before vs after the peak, and the interaction between the two. The direction of the effect, however, must be determined from the plots (below).

After looking at the plots, we see that when we account for the number of flowers, deviation has a greater negative effect on the late blooming Potentilla individuals.

4.3.2 Pollen Limitation Model

Below is the mixed effects model for the pollen limitation treatments.

#pot.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), family = "nbinom2", data = pot) #Potentilla model for negative binomial data with glmmadmb package
#summary(pot.treat.glmmtmb) #summary of model output

#pot.counts.treat.simulation<-simulateResiduals(fittedModel = pot.treat.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
#plot(pot.counts.treat.simulation) #plotting simulation
#testDispersion(pot.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(pot.counts.treat.simulation) #checking for zero inflation

We detected a significant effect of the pollen limitation treatment on the total number of developed seeds for Potentilla individuals, but the Potentilla pollen limitation data were zero-inflated. We again included the zero-inflation formula into our models.

pot.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), ziformula=~1, family="nbinom2",data = pot) #Potentilla model for negative binomial data with glmmadmb package
summary(pot.treat.glmmtmb) #summary of model output
##  Family: nbinom2  ( log )
## Formula:          total.seeds ~ treat/X.flowers + (1 | site/plot.treat)
## Zero inflation:               ~1
## Data: pot
## 
##      AIC      BIC   logLik deviance df.resid 
##    480.5    494.4   -232.2    464.5       34 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance  Std.Dev. 
##  plot.treat:site (Intercept) 3.090e-09 5.559e-05
##  site            (Intercept) 1.307e-02 1.143e-01
## Number of obs: 42, groups:  plot.treat:site, 6; site, 4
## 
## Overdispersion parameter for nbinom2 family (): 1.79 
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              3.874253   0.573392   6.757 1.41e-11 ***
## treatopen-hand           0.636866   0.720335   0.884   0.3766    
## treatopen:X.flowers      0.055490   0.032033   1.732   0.0832 .  
## treatopen-hand:X.flowers 0.009497   0.029244   0.325   0.7454    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.5791     0.6073  -4.247 2.17e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pot.counts.treat.simulation<-simulateResiduals(fittedModel = pot.treat.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
plot(pot.counts.treat.simulation) #plotting simulation

testDispersion(pot.counts.treat.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.78442, p-value = 0.312
## alternative hypothesis: two.sided
testZeroInflation(pot.counts.treat.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.97784, p-value = 1
## alternative hypothesis: two.sided

The data are no longer zero-inflated, and we still detected a marginally significant effect of pollen limitation in the Potentilla individuals. Potentilla individuals are slightly pollen limited.

4.3.3 Figures

Below are the plots for Potentilla individuals.

ggplot(pot.open, aes(x=pot.open$relative.position, y=pot.open$total.seeds, group=plot.treat)) +
  labs(title="Effect of Potentilla Blooming Time on Total Developed Seeds", y="Total Developed Seeds",x="Days Relative to Population Peak Bloom") +
  geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
  geom_smooth()

ggplot(pot.open, aes(x=pot.open$relative.position, y=pot.open$total.seeds, group=pot.open$early.late)) +
  labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Potentilla Blooming Time on Total Developed Seeds") +
  theme_classic() +
  geom_point(aes(color=early.late),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
  geom_smooth(method = "glm", aes(color=early.late))

ggplot(pot.open, aes(x=pot.open$deviation, y=pot.open$total.seeds/pot.open$X.flowers, group= pot.open$early.late)) +
  labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom",title="Effect of Potentilla Deviation from Peak on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
  geom_smooth(method="glm",aes(color=early.late)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))

ggplot(pot.open, aes(x=pot.open$number.conspecifics, y=pot.open$total.seeds/pot.open$X.flowers, group= pot.open$plot.treat)) +
  labs(y="Total Developed Seeds per Flower",x="Conspecific Density", title="Effect of Potentilla Conspecific Density on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
  geom_smooth(method="glm",aes(color=plot.treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))

Again we plotted the pollinator activity period.

ggplot(pot,aes(x=as.Date.factor(pot$midpoint), y=pot$total.seeds/pot$X.flowers, group= pot$treat)) +
  labs(y="Total Developed Seeds per Flower",x="Date", title="Pollinator Activity Period and Total Developed Seeds in Potentilla") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=treat)) +
  geom_smooth(aes(color=treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Open or Hand-pollinated") +
  scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  annotate("rect",xmin=as.Date.character("2019-07-05"), xmax=as.Date.character("2019-08-24"),ymin=0,ymax=Inf,alpha=0.1) +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d")

5 Mixed Effects Models with Proportion Data

Only the Mertensia and Delphinium data included proportion of developed seeds to total seeds. When counting seeds for the Potentilla individuals, undeveloped seeds were difficult to distinguish from other elements of the carpel and were not included. Because the data are proportions and we are using the cbind function, we assume a binomial distribution.

5.1 Mertensia

First, we had to calculate undeveloped seeds to use the cbind function for our proportion of developed seeds.

mert.binom<-fitdist(mert.open$prop.dev.seeds, "norm",start = NULL) #fitting Mertensia data to binomial distribution with normal errors using fitdistplus
plot(mert.binom) #plotting to see the fit

mert.open$undev.seeds<-(mert.open$total.seeds-mert.open$dev.seed) #creating an undeveloped seed column for cbind

mert$undev.seeds<-(mert$total.seeds-mert$dev.seed) #creating an undeveloped seed column for cbind

5.1.1 Main Model

Below is the mixed effects model for the proportion of developed seeds. The individuals are still only limited to the open treatment.

Mert.prop.glmmtmb<-glmmTMB(cbind(dev.seed,undev.seeds) ~ flowers:(deviation * early.late) + flowers:number.conspecifics/plot.treat + (1|site), family = binomial, data = mert.open) #Mertensia model with proportion of developed seeds using the glmmTMB package
summary(Mert.prop.glmmtmb)
##  Family: binomial  ( logit )
## Formula:          
## cbind(dev.seed, undev.seeds) ~ flowers:(deviation * early.late) +  
##     flowers:number.conspecifics/plot.treat + (1 | site)
## Data: mert.open
## 
##      AIC      BIC   logLik deviance df.resid 
##    306.8    319.4   -146.4    292.8       38 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.412    0.6419  
## Number of obs: 45, groups:  site, 4
## 
## Conditional model:
##                                               Estimate Std. Error z value
## (Intercept)                                  7.077e-01  3.599e-01   1.966
## flowers:deviation                           -1.950e-03  9.712e-04  -2.008
## flowers:early.latelate                      -5.714e-03  6.106e-03  -0.936
## flowers:number.conspecifics                  7.322e-05  7.290e-05   1.004
## flowers:deviation:early.latelate             3.983e-04  1.719e-03   0.232
## flowers:number.conspecifics:plot.treatmanip  4.777e-05  9.475e-05   0.504
##                                             Pr(>|z|)  
## (Intercept)                                   0.0492 *
## flowers:deviation                             0.0446 *
## flowers:early.latelate                        0.3494  
## flowers:number.conspecifics                   0.3152  
## flowers:deviation:early.latelate              0.8168  
## flowers:number.conspecifics:plot.treatmanip   0.6141  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mert.prop.simulation<-simulateResiduals(fittedModel = Mert.prop.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.prop.simulation) #plotting simulation

testDispersion(mert.prop.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0237, p-value = 0.856
## alternative hypothesis: two.sided
testZeroInflation(mert.prop.simulation)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.97847, p-value = 1
## alternative hypothesis: two.sided

We found a significant effect of deviation from the population peak on the proportion of developed seeds in Mertensia individuals. The Mertensia proportion data were also not overdispersed or zero-inflated. After looking at the plots below, we see that when flower number was taken into account, Mertensia individuals that deviated from the population peak had higher seed set.

5.1.2 Pollen Limitation Model

Next, we created the pollen limitation model for the Mertensia developed seeds. The model includes site and plot treatment as random effects like the models above. Now the proportion of developed seeds is the response variable.

Mert.prop.treat.glmmtmb<-glmmTMB(cbind(dev.seed,undev.seeds) ~ treat/flowers + (1|site/plot.treat), family = binomial, data = mert) #model with pollen lim treatment and proportion of developed seeds using the glmmTMB package
summary(Mert.prop.treat.glmmtmb) #summary of model results
##  Family: binomial  ( logit )
## Formula:          
## cbind(dev.seed, undev.seeds) ~ treat/flowers + (1 | site/plot.treat)
## Data: mert
## 
##      AIC      BIC   logLik deviance df.resid 
##    773.4    788.6   -380.7    761.4       87 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance Std.Dev.
##  plot.treat:site (Intercept) 0.1327   0.3643  
##  site            (Intercept) 0.0132   0.1149  
## Number of obs: 93, groups:  plot.treat:site, 8; site, 4
## 
## Conditional model:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             0.736049   0.202477   3.635 0.000278 ***
## treatopen-hand         -0.094499   0.194890  -0.485 0.627759    
## treatopen:flowers      -0.002515   0.002858  -0.880 0.378807    
## treatopen-hand:flowers -0.001949   0.002921  -0.667 0.504652    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mert.prop.treat.simulation<-simulateResiduals(fittedModel = Mert.prop.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.prop.treat.simulation) #plotting simulation

testDispersion(mert.prop.treat.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0213, p-value = 0.784
## alternative hypothesis: two.sided
testZeroInflation(mert.prop.treat.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99522, p-value = 1
## alternative hypothesis: two.sided

We didn't detected a significant effect of pollen supplementation on the proportions of developed seeds for Mertensia individuals. We again checked to see if the data were overdispersed or zero-inflated, and they were not. In contrast to the seed count results above, Mertensia individuals may not be pollen limited.

5.1.3 Figures

Below are the plots for Mertensia proportion data.

mert.open<-mert.open[!is.na(mert.open$prop.dev.seeds),] #removing NA values from proportion column

ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$prop.dev.seeds, group=plot.treat)) +
  labs(y="Proportion of Developed Seeds",x="Days Relative to Population Peak Bloom",title="Effect of Mertensia Blooming Time on Proportion of Developed Seeds") +
  geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
  geom_smooth()

ggplot(mert.open, aes(x=mert.open$deviation, y=mert.open$prop.dev.seeds/mert.open$flowers, group= mert.open$early.late)) +
  labs(y="Proportion of Developed Seeds per Flower",x="Days Relative to Population Peak Bloom",title="Effect of Mertensia Deviation from Peak on Proportion of Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
  geom_smooth(method="glm",aes(color=early.late)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))

Despite a marginally significant interaction effect, the effect of deviation on early vs. late individuals was not strong for Mertensia proportions of developed seeds.

ggplot(mert.open, aes(x=mert.open$number.conspecifics, y=mert.open$prop.dev.seeds/mert.open$flowers, group= mert.open$plot.treat)) +
  labs(y="Proportion of Developed Seeds per Flower",x="Conspecific Density",title="Effect of Mertensia Conspecific Density on Proportion of Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
  geom_smooth(method="glm",aes(color=plot.treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))

Finally, we added pollinator activity.

ggplot(mert, aes(x=as.Date.factor(mert$midpoint), y=mert$prop.dev.seeds/mert$flowers, group= mert$treat)) +
  labs(y="Proportion of Developed Seeds per Flower",x="Date", title="Pollinator Activity Period and Proportion of Developed Seeds in Mertensia") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=treat)) +
  geom_smooth(aes(color=treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Open or Hand-pollinated") +
  scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  annotate("rect",xmin=as.Date.character("2019-06-19"), xmax=as.Date.character("2019-07-03"),ymin=0,ymax=Inf,alpha=0.1) +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d")

5.2 Delphinium

We had to find the total developed seeds and undeveloped seeds before creating the model. To do this, we combined the seed counts of the first flower in bloom with the rest of the flowers on the plant.

delph.binom<-fitdist(delph.open$prop.dev.seeds, "norm") #fitting delphinium data to binomial distribution with normal errors using fitdistplus
plot(delph.binom) #plotting to see the fit

5.2.1 Main Model

Below is the mixed effects model for the proportion of developed seeds. The individuals are still only limited to the open treatment.

Delph.prop.glmmtmb<-glmmTMB(cbind(dev.seed.total,undev.seed.total) ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), family = betabinomial, data = delph.open) #modeling Delphinium proportion of developed seeds with glmmTMB package
summary(Delph.prop.glmmtmb)
##  Family: betabinomial  ( logit )
## Formula:          
## cbind(dev.seed.total, undev.seed.total) ~ X.flowers:(deviation *  
##     early.late) + X.flowers:number.conspecifics/plot.treat +      (1 | site)
## Data: delph.open
## 
##      AIC      BIC   logLik deviance df.resid 
##    334.5    348.2   -159.2    318.5       33 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 5.607e-09 7.488e-05
## Number of obs: 41, groups:  site, 4
## 
## Overdispersion parameter for betabinomial family (): 4.85 
## 
## Conditional model:
##                                                 Estimate Std. Error z value
## (Intercept)                                   -0.2054242  0.3484833  -0.590
## X.flowers:deviation                            0.0159357  0.0112946   1.411
## X.flowers:early.latelate                       0.1346731  0.1556862   0.865
## X.flowers:number.conspecifics                  0.0003517  0.0017993   0.196
## X.flowers:deviation:early.latelate            -0.0073549  0.0337245  -0.218
## X.flowers:number.conspecifics:plot.treatmanip -0.0014028  0.0017679  -0.793
##                                               Pr(>|z|)
## (Intercept)                                      0.556
## X.flowers:deviation                              0.158
## X.flowers:early.latelate                         0.387
## X.flowers:number.conspecifics                    0.845
## X.flowers:deviation:early.latelate               0.827
## X.flowers:number.conspecifics:plot.treatmanip    0.428
delph.prop.simulation<-simulateResiduals(fittedModel = Delph.prop.glmmtmb) #creating simulation for Delphinium model
plot(delph.prop.simulation) #plotting simulation

testDispersion(delph.prop.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0908, p-value = 0.416
## alternative hypothesis: two.sided
testZeroInflation(delph.prop.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 6.4935, p-value = 0.04
## alternative hypothesis: two.sided

The seed set of Delphinium individuals wasn't significantly affected by relative blooming time and conspecific density. Our model was overdispersed to binomial errors, but not to betabinomial errors. Relative blooming time and conspecific density do not have an effect on the proportion of developed seeds in Delphinium individuals.

5.2.2 Pollen Limitation Model

We created the mixed effects model for the pollen supplemented Delphinium individuals. Like the Mertensia model above, site and plot treatment are random effects, the interaction between pollen limitation treatment and number of flowers is the fixed effect, and proportion of developed seeds is the response variable.

Delph.prop.treat.glmmtmb<-glmmTMB(cbind(dev.seed.total,undev.seed.total) ~ treat/X.flowers + (1|site/plot.treat), family = betabinomial, data = delph) #modeling pollen lim and proportion of developed seeds with glmmTMB package
summary(Delph.prop.treat.glmmtmb) #model summary
##  Family: betabinomial  ( logit )
## Formula:          
## cbind(dev.seed.total, undev.seed.total) ~ treat/X.flowers + (1 |  
##     site/plot.treat)
## Data: delph
## 
##      AIC      BIC   logLik deviance df.resid 
##    660.3    677.3   -323.2    646.3       77 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance  Std.Dev. 
##  plot.treat:site (Intercept) 3.190e-09 5.648e-05
##  site            (Intercept) 2.482e-01 4.982e-01
## Number of obs: 84, groups:  plot.treat:site, 7; site, 4
## 
## Overdispersion parameter for betabinomial family (): 5.69 
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)  
## (Intercept)              -0.53669    0.49501  -1.084   0.2783  
## treatopen-hand            0.59997    0.57195   1.049   0.2942  
## treatopen:X.flowers       0.19255    0.08753   2.200   0.0278 *
## treatopen-hand:X.flowers  0.10651    0.07554   1.410   0.1585  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
delph.prop.treat.simulation<-simulateResiduals(fittedModel = Delph.prop.treat.glmmtmb) #creating simulation for Delphinium model
plot(delph.prop.treat.simulation) #plotting simulation

testDispersion(delph.prop.treat.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0232, p-value = 0.824
## alternative hypothesis: two.sided
testZeroInflation(delph.prop.treat.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 2.8571, p-value = 0.008
## alternative hypothesis: two.sided

Delphinium proportions of developed seeds were significantly affected by pollen supplementation, and the data were overdispersed to the binomial, but not to betabinomial errors.

As shown with total developed seeds, Delphinium individuals are pollen limited in this system.

5.2.3 Figures

Below are the plots for the proportion of developed seeds in Delphinium individuals.

delph.open<-delph.open[!is.na(delph.open$prop.dev.seeds),] #removing NA values in proportion column
ggplot(delph.open, aes(x=delph.open$relative.position, y=delph.open$prop.dev.seeds, group=plot.treat)) +
  labs(title="Effect of Delphinium Blooming Time on Proportion of Developed Seeds", y="Proportion of Developed Seeds",x="Days Relative to Population Peak Bloom") +
  geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
  geom_smooth()

ggplot(delph.open, aes(x=delph.open$deviation, y=(delph.open$prop.dev.seeds/delph.open$X.flowers), group= delph.open$early.late)) +
  labs(y="Proportion of Developed Seeds per Flower",x="Days Relative to Population Peak Bloom",title="Effect of Delphinium Deviation from Peak on Proportion of Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
  geom_smooth(method="glm",aes(color=early.late)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))

ggplot(delph.open, aes(x=delph.open$number.conspecifics, y=delph.open$prop.dev.seeds/delph.open$X.flowers, group= delph.open$plot.treat)) +
  labs(y="Proportion of Developed Seeds per Flower",x="Conspecific Density",title="Effect of Delphinium Conspecific Density on Proportion of Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
  geom_smooth(method="glm",aes(color=plot.treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue")) + xlim(0,110)

Finally, we added pollinator activity.

ggplot(delph, aes(x=as.Date.factor(delph$midpoint), y=delph$prop.dev.seeds/delph$X.flowers, group= delph$treat)) +
  labs(y="Proportion of Developed Seeds per Flower",x="Date", title="Pollinator Activity Period and Proportion of Developed Seeds in Delphinium") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=treat)) +
  geom_smooth(aes(color=treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Open or Hand-pollinated") +
  scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  annotate("rect",xmin=as.Date.character("2019-06-27"), xmax=as.Date.character("2019-07-18"),ymin=0,ymax=Inf,alpha=0.1) +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d")

ggplot(mert, aes(x=treat, y=dev.seed/flowers)) +
  labs(y="Total developed seeds", x="Pollen limitation treatment") +
  theme_classic() +
  geom_col()

ggplot(mert, aes(x=treat, y=prop.dev.seeds/flowers)) +
  labs(y="Proportion developed seeds", x="Pollen limitation treatment") +
  theme_classic() +
  geom_col()

ggplot(delph, aes(x=treat, y=total.seeds/X.flowers)) +
  labs(y="Total developed seeds", x="Pollen limitation treatment") +
  theme_classic() +
  geom_col()

ggplot(delph, aes(x=treat, y=prop.dev.seeds/X.flowers)) +
  labs(y="Proportion developed seeds", x="Pollen limitation treatment") +
  theme_classic() +
  geom_col()

ggplot(pot, aes(x=treat, y=total.seeds/X.flowers)) +
  labs(y="Total developed seeds", x="Pollen limitation treatment") +
  theme_classic() +
  geom_col()
## Warning: Removed 25 rows containing missing values (position_stack).

6 Session Info

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lmerTest_3.1-3     scales_1.1.1       glmmTMB_1.0.2.9000 DHARMa_0.3.3.0    
##  [5] lubridate_1.7.9.2  MuMIn_1.43.17      fitdistrplus_1.1-3 survival_3.2-7    
##  [9] MASS_7.3-53        RColorBrewer_1.1-2 lme4_1.1-26        Matrix_1.3-2      
## [13] forcats_0.5.0      stringr_1.4.0      dplyr_1.0.2        purrr_0.3.4       
## [17] readr_1.4.0        tidyr_1.1.2        tibble_3.0.4       ggplot2_3.3.3     
## [21] tidyverse_1.3.0    knitr_1.31        
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-151        fs_1.5.0            doParallel_1.0.16  
##  [4] httr_1.4.2          numDeriv_2016.8-1.1 tools_4.0.3        
##  [7] TMB_1.7.18          backports_1.2.1     R6_2.5.0           
## [10] DBI_1.1.0           mgcv_1.8-33         colorspace_2.0-0   
## [13] withr_2.3.0         tidyselect_1.1.0    compiler_4.0.3     
## [16] cli_2.2.0           rvest_0.3.6         xml2_1.3.2         
## [19] labeling_0.4.2      digest_0.6.27       minqa_1.2.4        
## [22] rmarkdown_2.6       pkgconfig_2.0.3     htmltools_0.5.1    
## [25] fastmap_1.0.1       dbplyr_2.0.0        highr_0.8          
## [28] rlang_0.4.10        readxl_1.3.1        rstudioapi_0.13    
## [31] shiny_1.6.0         farver_2.0.3        generics_0.1.0     
## [34] jsonlite_1.7.2      magrittr_2.0.1      Rcpp_1.0.5         
## [37] munsell_0.5.0       fansi_0.4.1         lifecycle_0.2.0    
## [40] stringi_1.5.3       yaml_2.2.1          plyr_1.8.6         
## [43] grid_4.0.3          promises_1.1.1      parallel_4.0.3     
## [46] qgam_1.3.2          crayon_1.3.4        lattice_0.20-41    
## [49] haven_2.3.1         splines_4.0.3       hms_1.0.0          
## [52] pillar_1.4.7        boot_1.3-25         codetools_0.2-18   
## [55] stats4_4.0.3        reprex_0.3.0        glue_1.4.2         
## [58] evaluate_0.14       modelr_0.1.8        vctrs_0.3.6        
## [61] nloptr_1.2.2.2      httpuv_1.5.5        foreach_1.5.1      
## [64] cellranger_1.1.0    gtable_0.3.0        assertthat_0.2.1   
## [67] xfun_0.20           mime_0.9            xtable_1.8-4       
## [70] broom_0.7.3         later_1.1.0.1       iterators_1.0.13   
## [73] statmod_1.4.35      ellipsis_0.3.1      gap_1.2.2